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Introduction 

Algorithm: A set of rules for solving a problem 
in a finite number of steps. 

Development: The progression to a more effec- 
tive state. 

The past decade has seen considerable activity 
in algorithm development for the Navier-Stokes 
equations. This has resulted in a wide variety of 
useful new techniques. It would appear, however, 
that there is plenty of room for further improve- 
ments. That is to say, we are far from exhausting 
all possible sets of rules for these problems and it 
is highly probable that some remaining ones will 
be more effective than those we have now. 

It is foolish and even counterproductive to an- 
ticipate or set milestones for the detailed develop- 
ment of basic or even applied research. The his- 
tory of science tells us that we can expect some- 
thing to happen in any major field if active minds 
capable of original thinking are allowed to pose 
challenging problems and seek elegant solutions. 

What we can do is look backwards and find 
what we are doing now in a given area of science 
that was not anticipated ten years ago. Some ex- 
amples of this type for the numerical solution of 
the Navier -Stokes equations form the body of this 
paper. These are divided into two parts, one de- 
voted to the incompressible Navier Stokes equa- 
tions, and the other to the compressible form. 
The discussion is far from being comprehensive, 
and. in fact, the examples for the incompressible 
case are strictly limited to experience at NASA 
Ames. 


1. Incompressible Navier-Stokes Equations 

In the middle and late 70s much attention was 
paid to the direct solution of homogeneous tur- 
bulent flows with periodic boundary conditions, 
see Rogallo (1981). The grids used at that time 
were M 6 and the storage capacity of available 
computers was the limiting factor in the spa- 
tial resolution. The natural method to use for 
the numerical approximation of the space deriva- 
tives was the classical spectral method composed 
of finite Fourier series, and the algorithm used 
for implementation was the fast Fourier trans- 
form. The time advance was fully explicit so 
that all of the time and space scales could be re- 
solved as accurately as possible. However, even 
with explicit schemes, time advance of a spectral 
method requires a minimum of two memory loca- 
tions for every dependent variable at every point 
in the mesh. The standard third and fourth or- 
der Runge Kut t a methods both take a minimum 
of three locations, so the options at that time 
were to use an existing second order time march 
method or use a coarser mesh and reduce the 
space accuracy. This motivated Dr. Alan Wray 
(Ames Research (’enter, unpublished) to turn to 
the Runge-Kutta technique and search for a sub- 
set of high order methods that require a minimum 
(two location) amount of storage capacity. He 
was successful and his third-order method, used 
to time march the nonlinear equation, 


du . 

h = F{u - i) 
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has the predictor-corrector form 

v — u n + cxAtF(u n , t n ) 
v = u ri + AAtF(u n . /. n ) 

1 / = m + /<A/F(i), i n + AAt) (1.2) 

« = u + BAiF(it , / n + /I A/) 
w n+ i = u + iAtF(u , f n + (a + B)At) 

The value of u n is initially provided and stored. 
The value of u is then calculated and also stored. 
Then the value of ii is formed and overstores u n 
which is no longer needed. The process continues 
through u and ii, requiring at any intermediate 
step only two memory locations per dependent 
variable per mesh point. Finally, «■„ + 1 overwrites 
w , ii is discarded, and the cycle is repeated. 

There are four equations for the five coefficients 
in eq (1.2), so we have a one-parameter family of 
low storage, third order Runge Kutta methods. 
The four equations are 

a + 0 + ^=1 
(a + Bh + AP = 1/2 
(a + B) 2 ~i + A 2 p= 1/3 
AB~i = 1/6 

One particular solution is given by 

o = l/4, A = 8/15 , 0 = 0, B = 5/12 ,7 = 3/4 

(1-4) 

This method is still being used to time march 
codes for homogeneous turbulent flows. It- is a 
good example of an algorithm advance adding a 
new capability to an old concept. 

A major advance in algorithms for wall bounded 
turbulent simulations occurred in the early 1980's. 
At that time Leonard and Wray (1982) extended 
the concepts being used to compute homogeneous 
turbulent flows, to compute wall bounded turbu- 
lent flows in relatively simple geometries. Let U 
be the velocity vector, p the pressure, and v the 
kinematic viscosity. One solves the vector equa- 
tion expressing conservation of momentum, 

U t + U VU = -Vp + vV 2 V (1.5) 

under the constraints of continuity in the domain 
and no slip at the walls: 

V -t/ = o in D , U = 0 on d D (1.6) 


In homogeneous flows harmonic basis functions 
are used and these automatically satisfy the pe- 
riodic boundary conditions. Furthermore, it was 
easy to make the solutions solenoidal (V * U — 0) 
so the pressure term could be eliminated. The 
idea advanced by Leonard and Wray was to build 
the constraints (1.6) into the basis functions of 
a generalized spectral method for wall bounded 
flows, so that the constraints are automatically 
and exactly satisfied with each time advance, and 
do not need to be further enforced at each step in 
conjunction with (1.5). The solution is then ex- 
pressed as a linear combination of global vector 
“basis functions 1 ’ that each satisfy (1.6). Due to 
the constraints one needs to carry only two de- 
grees of freedom per spectral mode while other 
methods usually carry four, the three velocity 
components and the pressure. Thus, less com- 
puter storage is needed to achieve the better res- 
olution. For more details and further discussion 
see the paper by Leonard and Wray. 

Where they can be formed (this can be dif- 
ficult since they are geometry dependent), the 
choice of the generalized spectral basis functions 
greatly improves the numerical treatment of the 
spatial aspect of the problem. However, to get 
adequate resolution near the walls, the time inte- 
gration tends to be stiff due to the eigenstructure 
of t he viscous terms. Because of this Dr. Philippe 
Spalart (Ames Research Center, unpublished) de- 
vised the use of a hybrid time marching scheme 
which is implicit for the (linear) viscous terms 
and explicit for the (nonlinear) convection terms. 
Again because of low memory requirements he 
had been using existing 2nd order methods for the 
time march. However, he has recently extended 
Wray's Runge-Kutta technique and developed a 
hybrid method which is third order in time ac- 
curacy and still has the minimum (two location) 
storage requirements. Thus if we have the vector 
relat ion 


u t = 7V(tt) + Lu (1.7) 


w'here L and N are matrix operators that are 
linear and nonlinear, respectively, the sequence 
can be made third order accurate with the proper 
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choice of the coefficients in 

u = u n + A t\L ■ (aiu n + 0\ii) + 7 , N n ] 
u = u + A t[L ■ (a 2 u + /? 2 «) + I 2 N + ft N n ] 

W n -)-i = U + A t\L ■ + 03 u n+l) + 73 iV + C 2 -/V] 

( 1 . 8 ) 

The treat ment of the N terms is equivalent to 
that used in Wray’s scheme. Only one solution 
for the coefficients is known. This is given by the 
conditions that 


7 , = 0.7208762469, 
73 = 0.5778221005, 
02 = 0.0929740417. 
C, = 0.4724519312, 


72 = 0.4001233399, 

0i = 0.3703996503, 

0 3 = 0.1818702938, 

C 2 = 0.2263697562 

(1.9a) 


at\+0i = 71 , ot 2 +02 — 72+Ci > a 3+03 — 73 +C 2 

(1.9b) 

Equations (1.9b) assure that the length of each 
time substep is the same for both L and TV. The 
numerical stability bounds for the model equation 


u t = i\u — vu 


( 1 . 10 ) 


where iXu represents N(u) and — vu represents 
L * w, are AA / < \/3 and vA l < 47, which were 
quite adequate for Spalart’s purposes. 


A concept that is probably as “new” as one 
can find is the total variation diminishing (TVD) 
theory extended to finite differencing schemes by 
Harten (1983). In this work, Harten extended 
ideas concerning total variation properties of scalar 
hyperbolic differential equations to discrete dif- 
ferencing schemes. This was an important step 
forward in determining the “ground rules” for de- 
signing good shock capturing methods, although 
it is not clear how religiously they need be fol- 
lowed. A complete review of this subject would 
be a formidable task by any measure. We chose 
not to do this, but rather to take some of the orig- 
inal underlying concepts involved and present a 
new perspective which hopefully will inspire new 
ideas. 

Flux- Vector / Flux-Difference Splitting 

In this section, we discuss two basic philoso- 
phies in the construction of upwind algorithms 
for systems of equations: flux-vector and flux- 
difference splitting. Each has proved to be a pow- 
erful technique in extending the upwind schemes 
for scalar equations to systems of equations. By 
the late 1970’s, the theory for scalar hyperbolic 
equations was well established and several up- 
wind schemes for these equations had appeared 
in the literature. The model nonlinear conserva- 
tion equation 


2. Compressible Navier-Stokes Algorithms 

The development of compressible Navier-Stokes 
algorithms has also seen moments of inspiration 
in the last decade. We have taken several steps 
forward in the general development of algorithms. 
Some of these steps are via new concepts while 
most are the result of applying old concepts in a 
new' setting. An example of a concept that was 
newly introduced to practical application in the 
held of fluid mechanics is the flux-vector splitting 
developed by Steger and Warming (1979). This 
opened up a wide new range of applications of 
“upwind” algorithms for the Euler and Navier- 
Stokes equations. Similar concepts have evolved 
since that time, most noticeably flux-difference 
splitting algorithms. Both of these methods have 
succeeded in removing much of the “fine tun- 
ing” of parameters w'hich had plagued many al- 
gorithms previous to this time. A brief review of 
this work is given. 


u t + {f{u)) x = 0 ( 2 . 1 ) 

had been analyzed extensively bv Lax (1973) and 
others as an initial- value problem, yielding a fairly 
complete description of the equation and its so- 
lution. For smooth regions of initial data, (2.1) 
can be represented for a small time interval by its 
quasi-Iinear form 

df 

u t + a{u)u x — 0 a(u) — — - 

du 

While at discontinuities, an integral form of (2.1) 
describes the solution behavior. The quasi-Iinear 
form has characteristic solutions for small time 
intervals of the form: u(x.t) — Uq(x — at ), i.e. 
the solution is constant along the characteristic 
lines. d d x t = a. Upwind methods (more prop- 
erly referred to as characteristic oriented meth- 
ods) use this information by determining the lo- 
cal propagation direction, sgn(a), and adapting 
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differencing stencils accordingly. One of t he sim- 
plest upwind schemes using this strategy is the 
Cole-Murman scheme. This scheme can be writ- 
ten for the discrete mesh, u™ — u{jAx,nAt), as 




*r + t = 5(/r + .+jr)-;i 5 ir + i(»" + .-«r) < 2 - 2 > 


where + i 


ini — li. if u u j + i 

a(uj) otherwise 


This produces the following simple (and more rec- 
ognizable) schemes for cases in which a is of uni- 
form sign: 


= if « >0 

»r' = if «<° 

Obviously, if higher order accuracy is needed, 
then a more elaborate scheme needs to be con- 
structed. But even for the simplest schemes (the 
Cole-Murman scheme for instance) one can ask 
the following question: what is the simplest and 
most natural way to extend the scalar upwind 
schemes to systems of equations ? For the Eu- 
ler equations, Steger and Warming (1979) and 
van Leer (1982) answered this question with flux- 
vector splitting while Roe (1981), Osher(1981) 
and others answered with a flux-difference split- 
ting technique. To illustrate these methods, we 
consider the 1-D Euler equations 

Q ( + (9 t E(Q) = 0 (2.3) 

Here Q is the vector of conserved variables for 
mass, momentum, and energy while E is the cor- 
responding flux vector. Whenever needed, we as- 
sume the ideal-gas law as an equation of state. 

The basic notion in flux-vector splitting is to 
split the flux vector into two parts 

E = E + + E~ 

The components, E _ and E + , are to be chosen 
such that they can be forward and backward dif- 
ferenced, respectively. This choice is based on 


the assumption that if the individual vectors can 
be forward and backward differenced in a stable 
fashion, i.e., if 

Q , + E 7 = 0 (stable for forward differencing) 

Q 1 + E+ — 0 (stable for backward differencing) 

(2.4a) 

t hen the same differencing can be used for the full 
equation, 

Q t + E++E;=0 (2.4b) 

This turns out to be the case, albeit some reduc- 
tion in stability characteristics may be encoun- 
tered. For the van Leer splitting described be- 
low with first order explicit upwinding, van Leer 
(1982) mentions that this amounts to a limit of 
GFL < 1 for (2.4a) and CFL < ^ for the full 
scheme, (2.4b). 

Steger and Warming constructed a general class 
of flux-vector splittings for the Euler equations by 
exploiting the fact that the flux vectors are homo- 
geneous of degree one in the conserved variables. 
Euler’s identity then gives 

dE , x 

E = ylQ with A = (2.5) 

()Q 

To construct the splittings, they first diagonalized 
A, 

'a, 

X~ ] AX = A = A 2 

A3 

then split the diagonal matrix into nonnegative 
and nonpositive parts, i.e. 

A = A + + A* (2.6) 

They define the new flux vectors by using (2.5) 

E ± = X\ ± X~'Q = A ± Q (2.7) 

The splittings based on (2.6) are obviously not 

unique; Steger and Warming suggest several dif- 
ferent splittings satisfying (2.6), of which proba- 
bly the most frequently used is the “ ± splitting 
defined by 

\± _ A F. I A I 

For this ± splitting they were able to determine 
that the resulting flux vectors had Jacobians, 





and , with all positive and negative eigen- 
values, respectively (at least for 1 < 7 < |). 
This is remarkable since nowhere in the develop- 
ment is any effort made to insure this. Unfortu- 
nately, this ± splitting leads to flux vectors that 
do not vary smoothly near sonic and stagnation 
points, even though the correct solution behaves 
smoothly there, and this produces “glitches” in 
the numerical solution. Several “fixes” have been 
proposed to remedy this, see Buning (1983). 

Van Leer (1982) provided an alternate flux- 
vector splitting, which he devised using special 
Mach number polynomials to construct fluxes that 
remain smooth near stagnation and sonic points. 
His construction technique is quite different from 
that of Steger and Warming, in particular it is 
purely a vector construction (neither the Jaco- 
bian matrix nor its diagonalizing transforms is 
directly used). A reasonable question to ask is 
whether van Leer's flux vectors have an equiva- 
lent Steger- Warming representation via similarity 
transforms of A as in (2.7). We find that this is 
so by redefining new and direct calculation. 
In the case of van Leer’s splitting they are given 

by 


A + = 


Mi 



A“ = A - A + 


with 

__ -((« - c) 2 - c 2 (7+ 1))(u + c) 2 
^ 1 4c 3 (7T 1) 

_ (u - C )(h~ + 2c )( w + c ) 2 

4c 3 (7+ 1 ) 

__ (( 7— 1 )u 2 + (1-37)uc + 2(27+1)c 2 )(u+c) 2 
/X3 4c 3 (7+ 1) 


In general, we find that these entries are of 
no particular uniform sign, (i.e. A + may have 
negative diagonal entries). This is not too sur- 
prising since the van Leer splitting only requires 
that the Jacobian matrices of the split fluxes, 
and , have nonnegative and nonposi- 

tive eigenvalues, respectively. For illustration, we 
chose the state: p — . 9 , u = .5 and c — 1 . 1 . 
In this case, the van Leer splitting gives: ~ 

.5097, — —.2885, /x 3 = 1.5098, while the eigen- 
values of are calculated to be 0.,.5795, and 


1.6918. Thus it appears that defining splittings 
from ( 2 . 6 ) is certainly not a necessary condition. 
We have, in fact, considered other schemes which 
satisfy ( 2 . 6 ) yet fail to have eigenvalues of their 
Jacobians with signs consistent with (2.6). This 
is certainly an avenue for future investigation. 

Flux-difference splitting has also provided a 
useful technique for extending scalar upwind al- 
gorithms to systems of equations. These methods 
use Riemann solvers to calculate the interaction 
of neighboring cells by the exact or approximate 
solution of Riemann’s initial- value problem. The 
simplest explicit schemes for solving the Euler 
equations take the generic structure: 


q; 


n +! _ Q n V» , — h 

ZL + A+l Zt 

At Ax 


= 0 


(2.9) 


where h, , i is the numerical flux at the cell inter- 
face between the grid points j and j + 1 . The role 
of the local Riemann solver is to determine the 
numerical flux at every cell interface by examin- 
ing the neighboring conditions. The best known 
approximate Riemann solvers are those of Roe 
(1981) and Osher and Solomon (1982). Roe's Rie- 
mann solver is particularly popular because of its 
simplicity. Roe considered the exact solution to 
the linearized form of ( 2 . 3 ), 


Q< + ^4(Q l .Q 7? )Qi = 0 (2.10) 


with constant left and right states specified as 
initial data, 


Q = 



x < 0 , t — 0 
x > 0 , t — 0 


Here x = 0 corresponds to the local cell inter- 
face and A is the approximate Jacobian, obtained 
from a mean value linearization satisfying 

E r - E 1 ' = /l(Q , %Q ff )(Q ff - Q l ) (2.11) 


Equation (2.10) can be diagonalized, decoupled, 
solved exactly, then recoupled. This amounts to 
solving three linear (scalar) convection problems 
with step functions as initial data and constant 
convection velocities u, u + c, and u — c. Since the 
exact solution for each scalar problem is merely 
the translation in x of the original step function, 
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this results in a “shocks only" approximate Rie- 
mann solver; expansion fans, shocks, and con- 
tact discontinuities are all modelled as discontinu- 
ities. Unfortunately, this allows expansion shocks 
to form as solutions which must be precluded by 
special means (see llarten (1983) for examples). 

From the local solution, the numerical flux at 
the cell interface can be calculated. If we con- 
struct A + and A~ as in (2.7) and (2.8), then the 
numerical flux can be written with reference to 
the left or right states as (details can be found in 
Roe (1981,1986)) 

h(Q L ,Q R ) =E t + A~ (Q l ,Q r )(Q r - Q l ) 

=E r - a + (q l ,q r ){Q r - Q l ) 

( 2 . 12 ) 

Taking the average and applying the local solu- 
tion everywhere on the discrete grid, we obtain 
the final form ( \A\ = A+ - A~ ) 

h J+ i=2(E, + 1 +E J )-^(Q,.Q J+1 )|(Q, + 1 -Q J ) 

( 2 - 13 ) 

If we again look at cases in which the local eigen- 
values are of uniform sign (supersonic flow), we 
obtain the following conventional schemes 

Q n+1 = Q n - ~ (E"-E"_ ,) if [u,u-c,u+c] >0 
Ax J J 

Q n+1 = Q n --^(E” +1 -E”) if [it, U—C, u+c] <0 

If we contrast this with the Cole-Murman scheme 
(2.2), which can also be viewed as using a “shocks 
only” scalar Riemann solver, we see that (2.13) is 
a successful extension of a scalar upwind scheme 
to systems. 

We conclude this section by remarking that 
we have limited our discussion to 1-D inviscid 
flow. This is not really as restrictive as one might 
guess. Remarkable success has been attained by 
applying these ideas “dimension by dimension” 
to the two and three-dimensional Navier-Stokes 
equations, see Chakravarthy and Osher (1985) for 
some excellent examples. 

Basics of TVD Schemes for Scalar Equa- 
tions 

In this section, we briefly mention the key ele- 
ments used in the development of the TVD con- 
cept. More details as well as proofs can be found 


in the literature. The basic notion is to consider 
solutions, u(z,/.), of the single nonlinear conser- 
vation equation 


«/ + (/(»'))* = o, = «(w) (2.14) 

au 

In this case, we make the usual assumption that 
the solutions of interest are entropy-satisfying weak 
solutions with convex flux functions. In the sim- 
plest case, to avoid boundary conditions, the ini- 
tial value problem is considered in which the solu- 
tion is specified along the r-axis, u(:z,0) = go{z), 
either in a periodic or compact supported fash- 
ion. This problem has been treated extensively by 
Lax (1973). The solution can be depicted in the 
x — t plane by a series of converging and diverging 
characteristic straight lines. From the solution 
of (2.14), Lax provides the following observation: 
the total increasing and decreasing variations of 
a differentiable solution between any pair of char- 
acteristics are conserved . This means that in the 
absence of shocks the exact solution of (2.14) con- 
serves the total variation of the initial data in 
time. 


I {t + fo) — I [l- o) 



du(x, t) 
dx 


dx 


(2.15) 


Moreover, in the presence of shock waves it can 
be shown that the total variation of the exact so- 
lution actually decreases in time ( i.e. /(f-Ho) < 
I{t 0 )). A simple heuristic argument for this de- 
crease would be to consider solution data with a 
shock present, u(x,t), and consider reconstruct- 
ing the solution data at a previous time u(x,t — 
Ai). But using characteristics, it becomes quickly 
obvious that this cannot be done uniquely; infor- 
mation (solution variation) has been irretrievably 
lost in the shock formation. An equally important 
result from Lax’s observation comes from consid- 
ering a monotonic solution between two noninter- 
secting characteristics: between pairs of charac- 
teristics, rnonotonic solutions remain monotonic , 
(i.e. no new extrema are created). 


Although the properties described previously 
are those of the differential equation (2.14) and 
its solution, Harten developed a TVD criterion 
for numerical schemes by considering the discrete 
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form of (2.15) on a mesh u" = u(jAx,nAt), 
u[ui,u 2 ,« 3 , ...]. The discrete total variation in 
this case is defined as 

-f CO 

rv(u) = ^|tt >+1 -tt,| (2.16) 

— oc 

with a corresponding TVD condition 

TV (u n+1 ) < TV (u n ) (2.17) 

It is not difficult to show that this TVD condition 
is sufficient for monotonic data with bounded to- 
tal variation to remain rnonotonic, (to prove this, 
assume a new extremum is introduced and com- 
pute the new total variation). Although we will 
only use (2.17) to investigate conditions for con- 
structing TVD schemes, equation (2.17), along 
with consistency of the scheme with the differen- 
tial equation and satisfaction of the entropy in- 
equality, is enough to guarantee convergence to 
the weak solution(s) (see Harten (1983)). 

Equation (2.17) now provides us with an ad- 
ditional measure which will allow us to rule out 
many existing schemes which do not diminish the 
solution variation. More importantly, as we will 
see in the next section, it will be used to derive 
algebraic criteria which we can use to construct 
new TVD schemes. 


Matrix Interpretation of TVD Criteria 

Sufficient conditions for constructing TVD al- 
gorithms have been developed first by Harten 
(1983) and later in a more general form by Osher 
and Chakravarthy (1984), and Jameson and Lax 
(1984). In this section we demonstrate general 
sufficient conditions for TVD schemes. In devel- 
oping the criteria for general explicit schemes, we 
independently followed a path similar to that of 
Jameson and Lax, although their claim of neces- 
sary and sufficient conditions is generally agreed 
to be in error (Harten (1988) notes that this is 
the danger of using their compact indicial nota- 
tion). In the development of implicit schemes we 
depart from their strategy altogether and avoid 
the introduction of expansive operators. More 
importantly, we avoid the use of indicial notation 
in favor of a more compact matrix- vector nota- 
tion whenever possible. As a result, the natu- 
ral simplicity of constructing TVD schemes be- 
comes evident, and we are able to give another 


(and perhaps clearer) interpretation of sufficient 
conditions given by the previous authors. 

An important step in the development of TVD 
schemes arises from the form chosen to express 
these schemes. We find it convenient to use a 
generalization of the form used by Osher and 
Chakravarthy. Since the objective is to obtain 
bounds on the variation of u, the conservative dif- 
ference schemes are put in a general form which 
uses an “apparent” (p + q + 2) explicit and (p f + 
q f + 2) implicit stencil of the solution, u. 


«; +, + 




n+l 


* = ~P' 


=«’+ E <?(<%+. A, + .+,»" 

(2.18) 

where A J + + i — u y Because C and 

D are typically nonlinear functions of u at grid 
points which could be outside the apparent sten- 
cils, it should be clear that (2.18) is far from 
being unique. This nonuniqueness provides a 
large amount of freedom in designing schemes 
and is essentially the distinguishing feature of 
various schemes appearing in the literature. Al- 
though the algebraic details of putting a particu- 
lar scheme into the form of (2.18) are important, 
we are only interested the general principles in- 
volved in the construction of TVD schemes and 
refer the reader to the literature for specific de- 
tails. 


We begin our analysis of TVD schemes by rewrit- 
ing the discrete total variation in terms of the 
forward difference matrix, 2?, (shown here for a 
periodic domain) 


/ - 1 1 0 0 
0 -110 
0 0 -11 


0 \ 
0 
0 


: 


0 00-11 

0 0 0 0 - 1 ) 


from which j|2?u||j = TV (u) and the TVD condi- 
tion can be written 


||Pu"+ 1 || 1 <11011-11, (2.19) 
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In these expressions we are using the nota- 
tion for the conventional L\ vector norm, UvHi — 
\ v j\- Using the forward difference matrix, eq. 
(2.18) can be written 


[/ + 0MD\ u n+1 =[/ — (!— 0)MD] u n (2.20) 

Here M and M are matrix operators which can 
be nonlinear functions of u. This equation, with 
the free parameter 0, represents various explicit 
and implicit forms of the evolution ol u hi time. 
(We chose this particular form so that if M = M , 
then the scheme would be a generalized form 
of Marten's “linearized” implicit TVD scheme, 
Marten (1984).) One can also construct a scheme 
representing the timewise evolution of the varia- 
tion, Du. To do so multiply (2.20) from the left 
by D and regroup terms. 


[/ + 0DM}Du n + 1 = [/ — (1 — 6)DM}Du r ' (2.21) 

Symbolically this can be expressed in terms of the 
matrix operators R and L as 

L Pu n+1 — R Du n or Du n+1 = L~ X R Du n 

( 2 . 22 ) 

with 

L = [/ + ODM], £ = [/-(!- 6)DM] 


choosing the column whose sum is largest. Fur- 
thermore, we have the usual matrix norm inequal- 
ity ||£~ 1 £|| 1 < ||£~ 1 ||i||£||i, so in the more gen- 
eral case, it is clear from (2.23) that it is sufficient 
to show that ||£ -1 ||i < 1 and ||£||i < 1 (L\ eon- 
tractive). As we will see, these simple estimates 
will be enough to obtain the TVD criteria of pre- 
vious investigators. 

First consider the explicit operator R and mul- 
tiply it from the left by the summation vector s = 
[1,1,1,...,]]. It is clear that sD = [0, 0, 0, ...,0], 
so that that R has columns that sum to exactly 
unity, regardless of the particular choice of M. 
Because the L\ norm of R is simply the maxi- 
mum of the sum of absolute values of elements in 
the columns of £, it is obvious that a sufficient 
and necessary condition for ||£||i < 1 is for R 
to be a nonnegative matrix , ( i.e. R > 0). Thus 
for explicit schemes [6 — 0) to be TVD, we have 
the general sufficient condition that R be a non- 
negative matrix. We illustrate that this leads to 
Marten’s criteria for explicit schemes by consider- 
ing his particular explicit form of (2.18), (in his 
notation) 


.«+! __ 


< + c; +j a, + .« m -c7_,a j . 


The operator R in this case (again assuming a 
periodic domain) has the following banded struc- 
ture 


Next we take the L\ vector norm of eq. (2.22) 
and apply the matrix- vector norm inequalities. 
Thus 


||Pu n+, ||i <||£- 1 *||,||Pu"||, (2.23) 

and we find a sufficient, condition for the scheme 
to be TVD is that ||£- , £||, < 1 

Note that for the extremely restrictive case in 
which £ and R are not functions of u, the ba- 
sic definition of a matrix norm would guaran- 
tee that ||£ -1 £||i < 1 is both a necessary and 
sufficient condition for the scheme to be TVD. 
(Many monotone schemes would be included in 
this class.) Recall that the L\ norm of a matrix 
is obtained by summing the absolute value of el- 
ements of individual columns of the matrix and 


je= 


o 

o 
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3 2 

0 

0 


0 

c + 

1 + 

1 - c\ , 

J+2 




- c 




0 

0 

°>+i 


0 

0 


(2.24) 


We need only require that this matrix be non- 
negative to immediately arrive at Marten’s crite- 
ria: 


t7 + , 

> 0 

J+ 2 


c~ , 

> 0 

J+2 
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dominant, that is, 


i ~ c i + i 


c 7+i 


For the general form of £, we can construct 
the matrix in the same fashion and arrive at the 
same conditions given by Harten, Jameson and 
Lax, and Osher and Chakravarthy by requiring 
that this resultant matrix be nonnegative. 

Perhaps the more interesting use of a matrix 
interpretation comes when considering implicit 
schemes. Sufficient conditions would be to show 
that both £ _1 and R are L\ contractive. We 
have shown sufficient conditions for constructing 
< 1 . We now consider conditions for mak- 
ing nr-'il. < 1 . From the previous development, 
one way to do this would be to show that £ 1 is 
a nonnegative matrix with columns that sum to 
unity. At that point the development would be 
the same as previously discussed. This turns out 
to be a simple task and using some well known 
results from matrix theory, we can determine al- 
gebraic sufficient conditions on £. 

Note that in the following discussion, we as- 
sume that £ is invertible, but after we have found 
a TVD criterion we will see that this must be so. 
First, we show that columns of £ _1 must sum to 
unity. We use the same trick of premultiplying £ 
by the summation vector, s = [1,1,1, ..., 1]. 

s£ = s — ► s — s£ _ 1 


Therefore the columns of £ _1 sum to unity. We 
need only find conditions on £ to make its in- 
verse nonnegative, but from matrix theory we 
know that a matrix whose inverse is nonnegative 
( £ _1 > 0 ) defines a monotone matrix. There- 
fore a sufficient condition would be that £ is a 
monotone matrix. This is not particularly use- 
ful in itself, but a well known theorem from ma- 
trix theory allows us to develop a TVD criterion. 
Sufficient conditions for £ monotone can be ob- 
tained from the theory for diagonally dominant 
M-matrices, a specific type of monotone matrix 
with positive diagonal entries and negative off- 
diagonal entries. To make this point clear we 
summarize a proof which appears in several books 
on matrix theory (see Lancaster, pp. 531-532 or 
Ortega, pp. 53-54). We begin by defining a real 
nxn matrix, a tJ , and assume that a %% > 0 for each 
i and a tJ < 0 whenever i / j. If A is diagonally 


a, t > ^ ^ M, 2 — 1,2 

then A is an M-matrix. To prove this, we first 
let D = diag[«ii , 033 , ..., a nn \ and define B — 

/ - D~ x A. Note that B has zero elements on the 
main diagonal and that B > 0. Also the fact that 
A is diagonally dominant implies that 

n 

^|6 tJ |<l, i = l,2,...,n. 
j = i 

It follows immediately from Gersgorin’s theorem 
that the maximum eigenvalue of B is less than 
one ( hb < 1 )• Now we have D~'A = I - B, 
and [7 - B ] -1 can be Neumann expanded into 

\I-B\- 1 =I + B + B 2 i B 3 + ... 

Since B > 0 , we conclude that [7 - B}~' > 0. 
It follows that 7> _ M is an M-matrix and conse- 
quently that A is an M-matrix. 


Therefore, sufficient conditions for £ to be mono- 
tone are that L be a diagonally dominant M- 
matrix, i.e. diagonally dominant with positive el- 
ements on the diagonal and negative off-diagonal 
elements. Also note that because of the diagonal 
dominance, we now can guarantee invertibility of 
£ as mentioned earlier. Again, we can recover 
the results of other investigators from these con- 
ditions. We illustrate this using Marten's implicit 
form 


In this case, £ takes the general structure 


£= 


-D~ 


1-2 

0 

0 


0 

'+ D Ui +D tn 


D 


J+2 


D U 


0 

0 
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To obtain Harten ’s TVI) criteria for the im- 
plicit scheme, we need only require that this be 
an M-matrix to obtain the following conditions, 
as did Harten 



> 0 

J + i 


ir , . 

> 0 

J + 2 



We conclude this section by noting the under- 
lying conceptual simplicity. Once the schemes are 
placed in the form of (2.21), then sufficient condi- 
tions become very simple and naturally give rise 
to the basic concepts of non negative and M - ma- 
trices. 

3. Concluding Remark 

Looking back over the last ten years, we can 
see that ten years ago it would have been correct 
to say: “There will be considerable advances in 
algorithm development in the next decade. ” We 
believe it is reasonably safe to make the same 
statement at this time. 
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